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INTRODUCTION 

With  the  use  of  adaptive  methods  to  solve  time-dependent  partial  differen¬ 
tial  equations,  there  exists  a  requirement  to  compute  solutions  on  moving  nonu¬ 
niform  grids.  There  is  also  a  requirement  to  estimate  the  local  discretization 
error  as  feedback  to  modify  or  refine  the  mesh.  In  this  report,  we  discuss  the 
MacCormack  finite  difference  scheme  and  a  Richardson  extrapolation-based  error 
estimation  procedure  that  was  used  in  the  adaptive  algorithm  of  Arney  (ref  1) 
and  Arney  and  Flaherty  (refs  2,3)  to  solve  time-dependent  hyperbolic  systems  in 
one-  and  two-space  dimensions.  Examples  of  other  adaptive  methods  with  these 
requirements  are  Rai  and  Anderson  (ref  4),  Adjerid  and  Flaherty  (ref  5),  Bell 
and  Shubin  (ref  6),  and  Davis  and  Flaherty  (ref  7). 

Finite  difference  methods  use  a  mapping  to  transform  the  time  and  space 
variables  from  a  moving  nonuniform  mesh  to  a  stationary  uniform  mesh.  The 
method  used  to  compute  the  metrics  of  this  transformation  must  be  carefully  cho¬ 
sen  in  order  to  preserve  the  stability,  conservation,  and  accuracy  of  the  scheme 
(cf.  Thomas  and  Lombard  (refs  8,9)  and  Hindman  (refs  10,11)). 

The  MacCormack  finite  difference  scheme  has  had  wide  use  in  solving 
Eulerian  conservation  laws  for  fluid  dynamics.  The  recent  use  of  artificial 
viscosity  to  make  this  scheme  total  variation  diminishing  (TVD)  makes  it  more 
attractive  as  a  general  solver  for  problems  with  discontinuities  (cf.  Davis  (ref 
(12)  and  Roe  (13)).  In  the  following  section,  we  discuss  the  MacCormack  scheme, 
our  implementation  of  the  differencing  of  the  metric  terms,  adaptive  selection 
of  the  time  step,  and  the  TVD  artificial  viscosity  of  Davis  (ref  12).  The 
Richardson's  extrapolation-based  error  estimation  method  produces  a  pointwise 
approximation  of  the  local  discretization  error  which  can  be  used  to  construct 
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several  global  measures  of  the  discretization  error.  Next,  we  discuss  our  error 
estimate  and  its  implementation  on  a  moving  mesh.  Then  we  present  computational 
results  of  solutions  of  hyperbolic  problems.  Computations  were  performed  in 
one-  and  two-dimensions  on  stationary  uniform  and  moving  nonuniform  grids. 
Lastly,  we  discuss  the  utility  of  our  methods,  the  computational  results,  and 
future  work. 

SOLUTION  SCHEME 

Consider  the  hyperbolic  vector  systems  of  conservation  laws  in  two-space 
dimensions 

-•  «♦  -#  — » 

ut  ♦  fx(x»y»u,t)  +  gy(x,y,u,t)  =  0  ,  (x,y)  e  D,  t  >  0  (1) 

u(x.y.O)  =  u0(x,y),  (x,y)  e  D  U  3D,  (2) 

with  appropriate  well-posed  conditions  on  the  boundary  3D  of  a  rectangular 
domain  D. 

We  chose  to  implement  the  MacCormack  finite  difference  scheme  for  hyper¬ 
bolic  problems  because  of  its  general  applicability.  The  MacCormack  scheme, 
like  most  higher-order  methods,  will  suffer  a  reduction  in  order  on  a  moving 
nonuniform  grid.  Despite  this  fact,  proper  mesh  moving  and  node  placement  by  an 
effective  adaptive  procedure  provide  enough  efficiency  and  accuracy  to  compen¬ 
sate  for  this  order  reduction. 

MacCormack  Scheme 

In  order  to  discretize  Eq.  (1),  we  introduce  a  transformation 

5  =  ?(x,y,t) ,  n  =  n(x,y,t),  t  =  t,  (3) 

from  the  physical  (x,y,t)  domain  to  a  computational  (£,/?, t)  domain  where  a  uni¬ 
form  rectangular  grid  will  be  used.  Under  this  transformation,  Eq.  (1)  becomes 
-*  -*  -•  -*  -•  -* 

UT  +  +  <Jnnt  +  f$*x  +  frjtlx  +  9^y  +  OpHy  =  0  (4) 
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The  transformation  metrics  (^x'^y'^t'Hx'^y'^t)  are  related  to  the  metrics 
(x^ , Xjj , xT #y^»yrj»yT)  by  the  identities 

.  .  yn  f  _  -*n  f  i'ii’irAi'iul  n  -  2i 

-  j"  '  5y  -  -j  /  ?t  -  3  '  ^X  _  j  ' 

xe  (yjXT-x.yT) 

ny  =  3-  .  nt  =  -----3 -  .  J  =  ynx{  -  Xny^  (5) 

Using  Eq.  (5)  in  Eq.  (4)  gives 

;T .  ;f  ias.jssr»  * >  .  ;t  ?  ♦  ?,  «^,  *  ;,<:?>  *  s  ?  - .  <*> 

This  equation  can  be  rewritten  in  another  form  in  the  original  transformation 

metrics  by  further  substitutions  of  Eq.  (5)  into  Eq.  (2)  as 
— *  -♦  — *  ^ 

UT  -  U£(xT£x+yT£y)  -  un(xTnx+yTny)  ♦  +  fnIIx  +  9£*y  +  Q^y  =  0  (7) 

Some  authors  (cf.  Hyman  (ref  14)  and  Thompson  (ref  15)  prefer  to  write  this 

equation  in  still  another  form  as 

-*  -*  -*■  -*■  -*  -* 

ut  +  +  ^t  +  f^x  +  9^y  +  9rj1y  *  0  <8) 

A  uniform  space-time  grid  having  mesh  spacing  A (  x  4i|  x  at  is  introduced 

onto  the  computational  domain.  The  finite  difference  solution  at  (iA£,  mArf, 

~*n  - 

nAT)  is  referred  to  as  Uj  m.  A  similar  notation  is  used  for  the  fluxes  f  and  g 

and  the  metrics  in  Eq.  (5).  The  two-step  MacCormack  scheme  (ref  16)  uses  first- 

order  forward  temporal  and  spatial  difference  approximations  in  the  predictor 

step,  and  first-order  backward  differences  in  the  corrector  step.  The  predicted 
— n+1 

solution  U|m  satisfies 

-n+1  -n  -n  -n  n  -n  -n  n 

Ml ,m  =  U|fn,  "  C(u!+l,m"ui,m)  ^t>i ,m  +  (fi  +  l,m-fi,m)  ^x)i,m 

-n  -n  n  -n  — n  n 

+  (9j+i,m-9|,mH^y)i,m3  -  C(uf  ,m+l-ui,m)  (>?t)f  ,m 

+  (^l,m+l“fI,m)(nx)l,m  +  (9l,m+l~9|fmHTty)*fni] 
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The  metrics  (Sx)jfin,  etc.  are  computed  by  forward  differences.  The  corrected 
-n+1 

solution  Ujm  satisfies 

-n+1  j  -n  -n+1  -n+1  -n+1  n+1  -n+1  -n+1  n+1 

ul,m  =  2  >u!,m  +  Ml.m  “  C (Mi ,m-yi-l , m) (£t ) 1 ,m  +  (f 1 ,m) ( ?x) i ,m 

-n+1  -n+1  n+1  -n+1  -n+1  n+1 

+  (9|,m“9i-l,m)  (^y)  2, ml  "  £rj  ^Ml,m*Mi,m-l)  (ht)l,m 

-n+1  -n+1  n+1  -n+1  -n+1  n+1  , 

+  (-! ,m“fi ,m-l ) (nx ) i ,m  +  (9lfm"2l,m-l) (hy)i(m] )  d°) 

-n+1 

with  metrics  computed  by  backward  differences.  The  notation  fj  m  denotes 
-n+1 

f (Ul(m) »  etc.  The  use  of  first  forward  and  backward  difference  approximations 
for  the  metrics  implies  that  the  transformation  from  the  computational  to  the 
physical  domain  is  piecewise  tri linear  in  space  and  time  for  the  predictor  and 
corrector  steps.  Such  low  order  difference  approximations  are  responsible  for 
reducing  the  orders  of  the  MacCormack  scheme.  A  smoother  transformation  and  the 
use  of  higher-order  difference  approximations  of  the  metrics  could  be  used  to 
maintain  secord  order  accuracy. 

It  was  shown  by  Hindman  (refs  10,11)  that  this  differencing  of  Eq.  (6)  pro¬ 
duces  consistent  approximations.  Therefore,  a  uniform  flow  solution  is  main¬ 
tained.  Other  conservative  forms  for  the  transformed  equations  were 
investigated  by  Hindman  (ref  10)  and  found  to  be  less  efficient  or  needing  spe¬ 
cial  differencing  of  the  metrics  for  computing  consistent  approximations. 

Equation  (4)  is  conservative  on  a  moving  mesh.  We  show  this  for  a  one¬ 
dimensional  scalar  conservation  law  by  investigating  the  Rankine-Hugoniot  jump 
conditions  across  a  shock  discontinuity.  Consider  a  conservation  law  in  the 
form 

fl  (/  u  dx)  +  f (u)|  =  0  (11) 

®  ^  -00  -00 


(12) 
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The  jump  conditions  for  a  discontinuity  at  x  =  s(t)  satisfy 

*  *  ¥-} 

[u] 

where  [q]  indicates  the  jump  in  q,  and  s  =  --  denotes  the  shock  velocity  (ref 


A  conservation  law  on  a  moving  mesh  produced  by  a  transformation  of 
variables  to  a  uniform  stationary  mesh  satisfies 


(§-  +  /  _  UXS<^  +  f(u)|  =  0 


Assuming  the  existence  of  a  shock  discontinuity,  £  ■  r(T)  gives 


rx^tu]  +  /  (ux^)Td^  +  /  <ux£)Td£  +  f(u)j  =  0 


Using  the  chain  rule  provides  an  integrable  form 

rxrfu]  +  /  (uxT-f )*d£  +  /  {uxT-f)fdf  +  f(u)|  *  0  (IE 

-ao  r  -oo 

Integration  of  this  equation  gives  jump  conditions  in  the  computational  domain 


rx?[u]  -  [f]  +  xT[u]  =  0  (16) 

Since  s(t)  and  r(T)  are  related  by 

•  • 

s  =  r  x^  +  xT  (17) 

the  appropriate  jump  condition  in  Eq.  (12)  is  recovered. 

Variable  Time  Step 

The  explicit  MacSormack  scheme  has  a  stability  restriction  that  limits  the 
time  step  allowed  for  a  given  spatial  mesh.  For  efficient  computation,  the  time 
step  should  be  adaptively  set  close  to  the  maximum  allowed  by  the  Courant, 
Friedrichs,  lewy  theorem  (ref  18) 


At  = 


0.8 


(18) 


2/2  max(i|/,u) 

The  computational  mesh  has  been  selected  to  have  spacing  A£  =  An  =  1  and  the 
constant  0.8  provides  a  20  percent  margin  of  safety.  The  quantities  <p  and  w  are 
the  spectral  radii  of  one-dimensional  conservation  laws  on  moving  meshes,  i.e., 

<|>  =  max[(Xi-xT){x  +  (Pi-yT)?y]  (19a) 

w  =  max[(XT-xT)nx  +  (P i-yT)ny]  (19b) 

where  and  p-j  are  eigenvalues  of  fy(u)  and  gj(u).  These  eigenvalues  and  the 
metrics  in  Eqs.  (19a)  and  (19b)  are  evaluated  at  the  beginning  of  each  time 
step. 

Artificial  Viscosity 

The  MacCormack  scheme,  being  a  second  order  accurate  centered  scheme,  pro¬ 
duces  spurious  oscillations  near  discontinuities.  In  order  to  eliminate  or 
reduce  these  oscillations,  artificial  viscosity  or  dissipation  is  added  to  the 
solution  to  diffuse  the  discontinuity.  The  viscosity  is  often  problem- 
dependent,  and  considerable  "fine  tuning"  is  usually  needed  to  balance  the 
effects  of  the  spurious  oscillations  and  diffusion  (ref  19). 

We  use  an  artificial  viscosity  model  due  to  Davis  (ref  12)  which  is  not 
problem-dependent  and  only  requires  knowledge  of  i|>  and  w.  This  artificial 
viscosity  model  is  designed  to  convert  the  MacCormack  scheme  into  a  total 
variation  diminishing  (TVO)  scheme  in  one-dimension.  A  scheme  is  TVD  if  the 
total  variation  of  the  solution  to  an  initial  value  problem  is  nonincreasing  in 
time.  Recent  research  efforts  have  resulted  in  the  development  of  other  second 
order  accurate  TVD  schemes  (cf.  Osher  and  Chakravarthy  (ref  20)  and  Warming  and 
Beam  (ref  21)). 
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The  artificial  viscosity  of  Davis  (ref  12)  is  based  on  a  flux  limiter  that 
does  not  depend  on  explicitly  determining  the  upwind  direction  and,  with  a 
recent  modification  by  Roe  (ref  13),  does  not  affect  the  region  of  stability  of 
the  MacCormack  scheme.  Because  the  MacCormack  scheme  also  does  not  determine 
the  upwing  direction,  the  combined  use  of  the  MacCormack  scheme  and  Davis'  arti¬ 
ficial  viscosity  is  computationally  simpler  to  perform  than  many  other  TVD 
schemes.  The  artificial  viscosity  terms  are  calculated  from  the  solution  data 
at  the  beginning  of  the  time  step.  For  two-dimensional  problems,  separate 
dissipative  terms  are  calculated  in  the  £  and  n  directions,  resepectively . 


ERROR  ESTIMATION 

Accurate  a  posteriori  error  estimation  is  an  integral  part  of  an  adaptive 
software  system.  Error  estimation  can  be  the  most  expensive  part  of  an  adaptive 
procedure  and  an  important  goal  is  to  find  accurate  and  inexpensive  ways  of 
estimating  the  discretization  error  (cf.  Babuska,  et  al.  (refs  22,23)).  The 
error  estimation  technique  is  dependent  on  many  factors,  including  the  type  of 
solver  used  in  the  algorithm,  the  type  of  error  to  be  determined,  and  the  norm 
in  which  the  error  estimate  is  to  be  measured.  It  is  most  desirable  to  have  a 
procedure  that  provides  pointwise  estimates  of  the  error  which  can  then  be  used 
to  find  estimates  in  several  local  and  global  norms. 

Mesh  nonuniformity  affects  the  accuracy  and  convergence  of  numerical 
schemes  and  error  estimation.  The  effects  of  the  mesh  on  the  solution  scheme 
have  been  studied  by  Ciment  (ref  24),  Fritts  (ref  25),  Hoffman  (ref  26),  Osher 
and  Sanders  (ref  27),  Sanders  (ref  28),  and  Mastin  (ref  29).  Error  analysis 
seems  to  be  more  natural  and  further  developed  for  finite  element  schemes,  espe¬ 
cially  for  elliptic  c  id  parabolic  problems  (cf.  Adjerid  and  Flaherty  (refs  5, 
30),  Zienkiewicz  et  al.  (refs  31,32),  and  Babuska  and  Rheinboldt  (refs  33,34)), 
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where  relatively  inexpensive  local  calculations  are  used  to  provide  accurate 
global  spatial  error  estimates.  More  study  needs  to  be  done  to  find  less  expen¬ 
sive  and  more  accurate  error  estimates  for  finite  difference  schemes  for  hyper¬ 
bolic  problems. 

We  calculate  the  local  temporal  and  spatial  portions  of  the  discretization 
error,  using  an  algorithm  based  on  Richardson's  extrapolation.  Flaherty  and 
Moore  (refs  35,36)  and  Berger  and  Oliger  (ref  37)  also  use  Richardson's  extrapo¬ 
lation  to  estimate  error  on  uniform  meshes  for  their  local  mesh  refinement 
algorithms. 

Richardson's  Extrapolation  Error  Estimation 

We  develop  the  error  estimation  for  the  second  order  MacCormack  scheme  for 

a  linear  scalar  problem  in  two  dimensions.  Separate  pointwise  estimates  at  a 

t 

general  spatial  node  i,  at  time  t,  for  the  local  temporal  error  En-(t)  and  local 
s 

spatial  error  E^(t)  are  obtained  with  two  different  extrapolation  procedures. 

Consider  a  uniform  mesh  with  spacing  Ax  x  Ay  and  time  step  At.  Let  the 
exact  solution  at  node  i  and  time  t  be  denoted  as  u-j(t),  the  numerical  solution 
by  the  MacCormack  scheme  at  the  same  point  and  time  as  Un- (t;Ax,Ay,At) ,  and  the 
MacCormack  finite  difference  operator  as  L(Ax,Ay,At),  i.e., 

U-j  (t+At;Ax,Ay,At)  =  L ( Ax , Ay , At ) U -j  ( t ; Ax , Ay , At )  (20) 

Assume  that  the  local  error  has  a  Taylor's  series  expansion  of  the  form 

ui(t)  -  Ui (t;Ax,Ay,At)  =  At[cjAt2  +  C2Ax*  +  C3Ay*  +  ...]  (21) 

where  the  constants  c^,  C2»  C3,...  are  independent  of  the  mesh  spacing. 

To  estimate  the  spatial  component  of  the  error,  we  calculate  a  solution  on 
a  mesh  of  double  spatial  size  (2Ax  x  2Ay)  with  the  same  time  step  (At).  The 
local  error  on  this  mesh  satisfies 

u.j  ( t+At )  -  Ui  (t+At;2Ax,2Ay,At)  =  Af^At*  +  4c2Ax*  +  4c3Ay*  +  ...]  (22) 
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Subtracting  Eq.  (22)  from  Eq.  (21)  and  neglecting  higher  order  terms,  we 
obtain  an  expression  for  the  leading  term  in  the  spatial  portion  of  the  local 


error  for  the  MacCormack  scheme  on  the  Ax  x  Ay  x  At  mesh  as 

s 

Ei (t+At) :=  At[C2Ax2  +  C3Ay2] 

=  |  [U-j  (t+At;2Ax,2Ay,At)  -  l).j  ( t+At ;  Ax,  Ay ,  At )  ]  (23) 

Similarly,  an  estimate  of  the  temporal  portion  of  the  local  error, 
t 

E^t+At),  can  be  calculated  by  computing  another  solution  on  the  Ax  x  Ay  spatial 

mesh  using  two  time  steps  of  At/2,  subtracting  this  result  from  Eq.  (21),  and 

retaining  the  leading  order  term  as 

t 

Ei (t+At) :=  At[c^At2] 

=  |  [Ui(t+At;Ax,Ay,2(|-)5  -  U-j (t+At;Ax,Ay,At) ]  (24) 

The  leading  term  of  the  local  error  at  node  i  and  time  t+At  is 

t  s 

E-j  (t+At)  *  Ej(t+At)  +  E^j  (t+At )  (25) 

There  are  several  disadvantages  to  this  technique  that  should  be  noted:  (1) 
the  error  cannot  be  calculated  for  nodes  on  or  adjacent  to  the  boundary;  (2)  the 
solution  must  be  smooth  enough  for  the  Cj,  03,  and  C3  to  exist;  (3)  the  error 
estimation  costs  approximately  three  times  more  to  compute  than  the  solution; 
and  (4)  the  mesh  must  be  uniform.  Equation  (25)  may  still  be  useful  as  a  mesh 
refinement  or  motion  indicator  even  in  situations  where  jumps  in  the  solution 
render  it  invalid  as  an  estimate  of  the  error. 

Richardson's  extrapolation  can  be  done  in  a  more  classic  manner  provided 
that  we  are  willing  to  forego  separate  spatial  and  temporal  error  estimates. 
Thus,  the  error  at  node  i  in  a  solution  on  a  mesh  having  spacing  Ax  x  At  is 
estimated  by  calculating  a  second  solution  on  a  mesh  with  spacing  Ax/2  using  two 
time  steps  of  At/2.  According  to  Eq.  (22),  the  local  error  on  this  mesh 
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satisfies 

u-j(t+At)  -  U, (t+At;Ax/2,2(At/2) )  *  At[cjAtV4  +  c2Ax*/4  +  ...]  (26) 

Subtracting  Eq.  (26)  from  Eq.  (22)  and  neglecting  higher  order  terms,  we 
can  obtain  error  estimates  for  either  U, (t+At;Ax,At)  or  U-j  (t+At;Ax/2,2(At/2) ) 
provided  that  node  i  is  common  to  both  meshes.  Our  adaptive  method  carries  the 
fine  grid  solution  forward  in  time;  thus,  we  estimate  its  error  as 

1 


Ei(t+At)  =  i  At(c1At2+c2Ax*) 


■  \  [U-j  ( t+At ;  5-  2  { f  -  )  )  -  ( t+At ;  Ax ,  At )  ] 


(27) 


Using  this  procedure,  the  error  can  now  be  calculated  at  nodes  adjacent  to 
boundaries.  Even  though  this  error  estimate  costs  four  times  more  to  compute 
than  the  solution,  we  only  incur  this  overhead  in  the  first  time  step.  No  addi¬ 
tional  cost  is  incurred  if  portions  of  the  mesh  have  to  be  refined  because  the 
solution  on  the  refined  mesh  has  already  been  computed  and  stored  while  esti¬ 
mating  the  error  for  the  coarser  parent  mesh. 

Error  Estimation  for  a  Moving  Nonuniform  Mesh 

Nonuniformity  of  the  mesh  changes  the  discretization  error  of  the 
MacCormack  scheme.  For  simplicity,  we  will  determine  this  error  and  analyze  its 
effects  on  the  Richardson  extrapolation  error  estimation  using  a  linear  scalar 
problem  in  one-space  dimension 


ut  +  buv  -  0 


(28) 


The  local  error  for  the  MacCormack  method  on  a  one-dimensional  moving  nonu¬ 
niform  mesh  is 

b  n+1  n 

U-j  ( t+At )  -  U-j  (t+At;Ax,At)  =  At[-  ^  (Axr  -Axj)uxx 


Axn  +  1 

-  Atb2(l  -  (-----)uxx)  +  ciAt2  +  c2Ax* ] 


(29) 
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where  Axj  and  Axr  are  the  mesh  sizes  on  the  left  and  right  node  i  at  time  step 

n+1  n+1 

n,  respectively,  and  Ax  *  max(Axr  ,Axj  ).  On  the  moving  nonuniform  mesh,  both 
the  temporal  and  spatial  error  components  contain  second  order  terms,  whereas 
the  error  on  a  uniform  mesh  is  third  order.  The  previous  analysis  can  be  used 
to  show  that  the  leading  component  of  the  temporal  error  is 


t  Axn+1 

Ej  (t+At) At[-  Atb*(l  -  ---r-JUxxJ 


*  2[Ui(t+At;AX,2{g-))  -  U,  ( t+At ; Ax , At )  ]  (30) 

Calculation  of  the  spatial  portion  of  the  error  is  more  difficult  since  the 
temporal  portion  of  the  error  does  not  cancel  upon  subtraction  of  solutions 
calculated  on  two  spatially  different  meshes.  We  overcome  this  difficulty  and 
also  greatly  simplify  the  procedure  in  two  dimensions  by  constraining  the  mesh 
to  maintain  double  size  increments  for  special  nodes  of  the  moving  coarse  mesh. 
This  constrained  grid  structure  consists  of  a  coarse  mesh,  shown  with  darker 
lines  in  Figure  1,  containing  properly  nested  fine  cells  created  by  binary  divi¬ 
sion  of  the  sides  of  the  coarse  cells,  shown  by  lighter  lines  in  Figure  1.  The 
vertices  of  the  coarse  cells  are  denoted  as  "independent  moving  nodes."  Error 
estimates  are  calculated  for  these  nodes.  The  remaining  nodes  in  the  mesh  of 
Figure  1  are  "dependent  moving  nodes"  which  must  be  moved  to  maintain  the 
constrained  grid  structure.  A  solution  is  computed  for  these  "dependent  moving 
nodes,"  but  no  error  estimate  is  obtained. 

For  the  "independent  moving  nodes,"  the  spatial  error  calculation  can 
proceed  as  for  a  uniform  mesh;  therefore,  the  local  spatial  error  estimate  is 

s  b  n+1  n 

^i (t+At)  *  At [-  j(Axr  -AXf )uxx] 

*  U-j  (t+At;  Ax,  At)  -  U-j  (t+At ;  2Ax ,  At ) 


(31) 
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Figure  1.  Spatial  structure  of  the  moving  coarse  mesh  (bold  lines)  with 
embedded  fine  mesh  (fine  lines)  used  for  the  error  estimation. 

The  above  analysis  extends  directly  to  two  dimensions;  hence,  we  have  a 
Richardson  extrapolation-based  procedure  of  estimating  error  on  a  moving  nonuni¬ 
form  grid.  In  practice,  we  test  the  need  for  local  uniformity  and,  if  found, 
use  formulas  in  Eqs.  (23)  through  (25)  to  compute  error  estimates. 

Error  estimation  for  systems  of  equations  involves  the  use  of  a  vector  norm 
at  node  i  and  time  (t).  The  examples  in  the  following  section  use  the  maximum 
norm,  i.e., 

E-j  (t)  :=  max  |  E-j ,  j(t)|  (32) 

1  <  j  <  N 

where  N  is  the  number  of  equations  in  the  system  and  E^j(t)  is  the  local  error 
estimate  for  the  jth  component  of  the  solution  vector  at  node  i. 

COMPUTATIONAL  EXAMPLES 

The  solution  and  local  error  estimation  procedures  are  applied  to  four 
examples.  In  Example  1,  we  demonstrate  the  capability  of  the  MacCormack  scheme 
with  Davis'  TVD  artificial  viscosity  on  a  moving  nonuniform  mesh.  In  Example  2, 
we  investigate  a  one-dimensional  problem  using  a  modified  form  of  the  error 
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estimate  in  Eqs.  (27)  and  (28).  Examples  3  and  4  illustrate  the  performance  of 
the  error  estimation  procedure  on  a  problem  having  a  smooth  solution  and  one 
with  a  jump  in  the  first  derivative,  respectively.  We  investigate  the  accuracy 
and  convergence  of  the  local  error  estimator  by  determining  an  effectivity  index 


9  = 


at  a  fixed  time  t  for  several  different  meshes  and  different  adaptive  strate¬ 
gies.  Here  e  and  E  are  the  exact  and  estimated  errors,  respectively.  The  Lj 


norm. 


IIEIIi  :=  //  Edxdy 

is  obtained  by  assuming  E  to  be  a  piecewise  constant  function. 

Example  1 

Consider  the  initial-boundary  value  problem 

ut  -  yux  +  xuy  =  0,  t  >  0,  -1.2  $  X  «  1.2,  -1.2  (  y  U.2 


(33b) 


u(x,y,0)  =  - 


,  if  (x  -  *)*  ♦  1 . 5y2  »  J- 


1  -  16 ( (x  -  ±)  +  1.5y*)  ,  otherwise 


u( 1 . 2,y, t )  =  u(-1.2,y,t)  =  u(x,-1.2,t)  =  u(x,1.2,t)  =  0 


u(x,y,t)  =  - 


0  ,  if  C  <  0 


C  ,  if  C  >  0 


where 


C  =  1  -  16((xcost  +  ysint  -  -)  +  1.5(ycost  -  xsint)2)  (38) 

Equations  (37)  and  (38)  represent  a  moving  elliptical  cone  rotating  coun¬ 
terclockwise  around  the  origin  with  period  2n.  This  problem  was  proposed  as  a 
test  problem  by  Gottlieb  and  Orszag  (ref  38)  and  was  used  as  a  test  problem  in  a 
survey  by  McRae  et  al.  (ref  39). 


row* 


We  show  the  sequence  of  meshes  that  were  generated  at  t  *  0,  1.6,  and  3.2 
using  the  adaptive  mesh  moving  method  of  Arney  and  Flaherty  (ref  2)  in  Figures 
2,  3,  and  4,  respectively.  Arney  and  Flaherty's  mesh  moving  method  (ref  2)  uti- 
lizes  the  error  estimates  (see  Error  Estimation  section  of  this  report)  to  con¬ 
centrate  the  mesh  in  the  high-error  region  beneath  the  cone  and  to  follow  it  as 
it  rotates.  It  also  increases  the  accuracy  of  the  solution  and  reduces  oscilla¬ 
tions  in  the  wake  following  the  cone.  However,  small  oscillations  are  still 
present.  Next,  we  solve  this  problem  with  the  same  moving  mesh  technique  by 
using  Davis'  artificial  viscosity  (ref  12)  with  the  MacCormack  scheme.  Surface 
and  contour  plots  of  solutions  with  and  without  artificial  viscosity  are  shown 
in  Figures  5  and  6.  There  is  no  artificial  wake  behind  the  cone  when  artificial 
viscosity  is  used.  However,  the  artificial  viscosity  slightly  diffuses  the 
cone,  widening  its  base  and  reducing  its  peak  from  1.0  to  0.88. 
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Figure  2.  Initial  mesh  for  Example  1. 


Figure  5.  Contour  plots  of  the  solutions  of  Example  1  on  a  moving  mesh 
without  artificial  viscosity  (left)  and  with  artificial 
viscosity  (right)  at  t  *  3.2. 


MS 


kr»> 


wiiiiMi 


Figure  6.  Surface  plots  of  the  solutions  of  Example  1  on  a  moving  mesh 
without  artificial  viscosity  (top)  and  with  artificial 
viscosity  (bottom)  at  t  »  3.2. 
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Example  2 

We  consider  an  application  of  the  direct  Richardson's  extrapolation  error 
estimation  procedure,  Eq.  (28),  to  the  one-dimensional  linear  scalar  equation 

ut  +  ux  =  0,  t  >  0,  0  <  x  <  0.8  (39) 

with  initial  and  Oirichlet  boundary  conditions  specified  so  that  the  exact  solu¬ 
tion  is 

u(x,t)  =  |  [1  -  tan  h  100  (x-t-0.2)]  (40) 

This  solution  is  a  relatively  steep  wave  that  moves  at  unit  speed  across  the 
domain. 

We  solved  this  problem  for  one  time  step  on  seven  different  uniform  meshes 
having  N  computational  cells  per  time  step  in  order  to  investigate  accuracy  and 
convergence  of  the  error  estimate.  Table  I  shows  the  results  obtained  from 
these  calculations. 

TABLE  I.  EXACT  ANO  ESTIMATED  ERRORS  FOR  DIFFERENT  MESH  SIZES 
FOR  EXAMPLE  2 


At 

N 

Exact  Error 

Hell  ^ 

Estimated  Error 

II E II  i 

Effectivity  Ratio 

9 

0.1 

8 

0.447  x  10~z 

0.351  x  10"’ 

7.60 

0.05 

16 

0.234  x  10-* 

0.132  x  10-1 

5.55 

0.025 

32 

0.106  x  10"z 

0.236  x  10-8 

2.22 

0.0125 

64 

0.257  x  10'3 

0.138  x  10“3 

0.54 

0.00625 

128 

0.294  x  10"« 

0.380  x  10"* 

1.29 

0.00312 

256 

0.303  x  10"* 

0.453  x  10-5 

1.49 

0.00156 

512 

0.538  x  10-* 

0.661  x  10-» 

1.23 

m 


We  also  solved  this  problem  using  Arney  and  Flaherty's  adaptive  local 
refinement  procedure  (ref  3)  on  a  base  mesh  having  Ax  *  At  *  0.1  with  a  local 
error  tolerance  of  1/128.  The  mesh  created  by  the  local  refinement  algorithm 
is  shown  in  Figure  7  and  the  solutions  computed  at  each  base  time  step  are  shown 
in  Figure  8. 


Figure  7.  Adaptive  mesh  of  Example  2. 


Figure  8.  Solutions  at  base  time  steps  for  the  adaptive  local 
refinement  procedure  of  Example  2. 

The  adaptive  composite  mesh  of  Figure  7  shows  a  distinct  pattern  associated 
with  using  the  MacCormack  scheme  with  our  local  refinement  strategy.  Spacious 
oscillations  of  the  solution  on  the  base  mesh  cause  several  levels  of  refinement 
which  drastically  reduce  the  base  mesh  spacing  at  the  beginning  of  each  base 
time  step.  However,  once  these  oscillations  have  been  controlled,  the  need  for 
refinement  is  reduced  at  the  later  stages  of  the  adaptive  procedure.  This 
situation  could  be  alleviated  by  including  an  artificial  viscosity  model  with 
the  MacCormack  scheme. 
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Example  3 

Consider  the  linear  scalar  hyperbolic  differential  equation 

ut  +  2ux  +  2Uy  =  0,  t>0,  0.2  $  x  $  1.2,  0  >  y  $  1 

with  initial  conditions 

(1-tan  h  3 ( x-0 . ly+0 . 1 ) ) 
u(x,y,0)  =  - - "2 - L - -- 


(41) 


:42: 


and  with  Dirichlet  boundary  conditions  specified  so  that  the  exact  solution  of 
this  problem  is 

M-tan  h  Vv-n  1\/-1  fl-t-+n  1\\ 

(43) 


u (x, y , t )  =  il:S3D-b-2i5:2iIy:li§J+9iil2 


This  solution  is  a  smooth  wave  that  moves  at  an  angle  of  45  degrees  across 
the  domain.  The  problem  was  selected  to  show  the  convergence  and  accuracy  of 
Richardson's  extrapolation  error  estimation  procedure  (Eqs.  (28)  and  (29)).  We 
solve  Eqs.  (41)  and  (42)  for  one  time  step.  At  *  0.012,  on  eight  different 
meshes.  The  mesh  strategy  of  each  calculation  is  described  as  follows: 

1.  a  stationary  uniform  (10  x  10)  rectangular  mesh 

2.  a  stationary  uniform  (20  x  20)  rectangular  mesh 

3.  a  stationary  uniform  (40  x  40)  rectangular  mesh 

4.  a  stationary  uniform  (60  x  60)  rectangular  mesh 

5.  a  stationary  (40  x  40)  mesh  of  nonuniform  quadrilateral  cells 

6.  a  moving  (20  x  20)  mesh  with  uniform  rectangles 

7.  a  moving  (20  x  20)  mesh  of  nonuniform  quadrilateral  cells 

8.  a  moving  (40  x  40)  mesh  of  nonuniform  quadrilateral  cells 

Table  II  shows  the  results  from  these  calculations  by  comparing  the  exact  errors 
and  the  effectivity  indices  for  the  eight  strategies. 
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TABLE  II.  EXACT  AND  ESTIMATED  ERRORS  FOR  DIFFERENT  MESH 
STRATEGIES  FOR  EXAMPLE  3 


Mesh  Strategy 
(From  Above) 


Exact  Error 
lie  14 

Estimated  Error 

II E II 1 

Effectivity  Ratio 

9 

0.0111 

0.0071 

0.64 

0.00370 

0.00318 

0.86 

0.000942 

0.000908 

0.96 

0.000367 

0.000368 

1.00 

0.000399 

0.000418 

1.04 

0.00136 

0.00124 

0.91 

0.000411 

0.000370 

0.90 

0.000167 

0.000156 

0.94 

Strategies  1  through  4  show  the  convergence  of  the  error  estimates  on  uni¬ 
form  meshes  as  the  number  of  nodes  increases.  These  errors  show  a  rate  of  con¬ 
vergence  of  0(Ax2,Ay2)  which  is  predicted  in  Eq.  (21).  Comparison  cf  the  errors 
of  strategies  3  and  5  shows  the  error  is  cut  in  half  by  computing  with  a  better 
nonuniform  stationary  mesh.  Further  comparison  of  strategies  5  and  7  shows 
another  reduction  of  error  by  half  when  the  mesh  is  properly  moved.  The  nonuni¬ 
formity  of  the  mesh  in  strategies  5,  7,  and  8  produces  little  change  in  the 
effectivity  of  the  error  estimation.  These  nonuniform  mesh  computations  indi¬ 
cate  a  convergence  rate  0(Axl-32(^yl.32) . 

Example  4 

Consider  the  linear  scalar  hyperbolic  differential  equation 

ut  +  ux  +  0.25uy  =0,  t  >  0,  0.2  $  x  $  1.2,  0  $  y  $  1  (44) 
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with  initial  conditions 


The  solution  of  this  problem  is  an  oblique  ramplike  wave  front  that  moves 
at  an  angle  of  14  degrees  across  the  domain.  The  solution  has  a  jump  in  its 
first  partial  derivatives  at  the  top  and  bottom  edges  of  the  wave  front.  We 
expect  some  difficulty  in  estimating  the  error  near  locations  where  the  deriva¬ 
tives  jump.  In  the  region  of  the  front  itself,  the  gradient  of  the  solution  is 
constant  and  there  is  no  error  in  the  solution  or  in  the  error  estimate. 

We  solved  this  problem  for  one  time  step,  At  =  0.015,  for  the  following  six 
mesh  strategies: 

1.  a  stationary  uniform  (12  x  12)  rectangular  mesh 

2.  a  stationary  uniform  (24  x  24)  rectangular  mesh 

3.  a  stationary  uniform  (48  x  48)  rectangular  mesh 

4.  a  stationary  uniform  (64  x  64)  rectangular  mesh 

5.  a  stationary  (24  x  24)  mesh  of  nonuniform  quadri lateral  cells 

6.  a  moving  (24  x  24)  mesh  of  nonuniform  quadrilateral  cells. 

Table  III  shows  the  results  of  these  strategies. 
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TABLE  III.  EXACT  AND  ESTIMATED  ERRORS  FOR  DIFFERENT  MESH  STRATEGIES  FOR 
EXAMPLE  4.  THE  ERROR  ESTIMATE  IS  ACCURATE,  BUT  THE  SOLUTION 
APPEARS  TO  BE  CONVERGING 


Mesh  Strategy 

Exact  Error 
Belli 

Estimated  Error 

HE  II 2. 

Effectivity  Ratio 

0 

1 

0.0058 

0.0016 

0.28 

2 

0.00275 

0.00110 

0.40 

3 

0.000866 

0.000479 

0.55 

4 

0.000400 

0.000222 

0.56 

5 

0.00144 

0.00078 

0.54 

6 

0.000720 

0.000349 

0.49 

The  results  are  once  again  as  expected.  The  error  estimate  of  this  problem 
with  a  jump  in  the  derivative  is  not  as  accurate  as  the  smooth  solution  of 
Example  3.  However,  the  error  estimate  still  shows  signs  of  converging  to  the 
exact  error  in  lj  for  the  uniform  meshes  of  strategies  1  through  4.  Once  again, 
the  better  nodal  placement  of  the  initial  mesh  by  the  mesh  generator  of  Arney 
(ref  1)  reduces  the  error  by  half  from  a  uniform  mesh.  Also,  moving  the  mesh  by 
the  method  of  Arney  and  Flaherty  (ref  2)  reduces  the  error  by  half  again. 

CONCLUSION 

We  have  shown  that  MacCormack's  finite  difference  scheme  and  error  estima¬ 
tion  based  on  Richardson's  extrapolation  can  be  used  on  moving  grids  with  local 
refinement.  With  proper  computation  of  the  transformation  metrics  and  the  use 


of  TVO  artificial  viscosity,  the  MacCormack  scheme  is  stable  and  is  able  to 
solve  problems  with  sharp  discontinuities. 


The  examples  we  have  presented  demonstrate  the  utility  of  these  methods  and 


also  point  out  their  shortcomings.  Of  particular  concern  is  the  lack  of  any 
error  estimation  near  the  boundaries,  the  poor  error  estimation  near  discon¬ 
tinuities,  and  the  need  to  constrain  the  mesh  to  obtain  any  accurate  error  esti 
mation.  These  problems  must  be  solved  in  order  to  effectively  utilize  this 
solution  scheme  and  error  estimation  procedure  with  an  adaptive  technique. 
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